#include <iostream>
#include <pqxx/pqxx> 
#include <fstream>
#include <sstream>
#include <math.h>
#include <ctime>

#define SIZE 46
#define YEAR 2011

using namespace std;
using namespace pqxx;



double get_sum_diff2( double[], double[]);
double get_sum_product(double[]);
double get_sum_product(double[], double[]);
double get_sum_product_diff_avg(double[]);
double get_sum_product_diff_avg(double[], double[]);

double get_Minimum_squared_Euclidean_distance( double[], double[] );
double get_Cosine_similarity( double[], double[] );
double get_Correlation_coefficient( double[], double[] );

void get_px_NDVI_serie( int, int, int, double, double[], int* );
void clear_array(double[], int);
void print_array(double[], int);

/* Make sure your input is is formated and named respectively

"groundtruth_2011.txt"
   . (int) number of ground truth pixel.
   For next (#) line:
   . (int) (int) line and sample of that pixel

"threshold.txt"
   . (int)       number of cloudy day allow to process
   . (int)       for cloud band03_surf_refl (by xiao's algorithm)
   . (double)    for MSED
   . (double)    for Cosine similarity
   . (double)    for Correlation coefficient

"range_to_process.txt"
   . (int) (int)     top-left (line sample)
   . (int) (int)     bottom-right (line sample)

*/

/*
Expected output is formated as below and named 

"log.txt"
   .(int) (int) (double) (bool) (double) (bool) (double) (bool)
   ...
      pixel_line pixel_sample MSED 1/0 CS 1/0 CC 1/0
   ...
   .(int) total pixels in range
   .(int) total pixels are removed by threshold for cloud

Output is included 3 files:
   - msed.txt
   - cs.txt 
   - cc.txt is mapped as padding rice field;
   example: 
      .(int)    number of paddy (pixel)
      .01001011100100
      .10010010100001
      .11010101101100
      .10110011100110

*/

int main(int argc, char* argv[])
{
   time_t show_start = time(0);
   // get groundtruth data position
   ifstream groundtruth;
   groundtruth.open("groundtruth_2011.txt");
      int number_of_gt_pixel;
      groundtruth >> number_of_gt_pixel;
      double gt_coordinate[2][number_of_gt_pixel];
      for (int i = 0; i < number_of_gt_pixel; ++i)
      {
         groundtruth >> gt_coordinate[0][i]; // line
         groundtruth >> gt_coordinate[1][i]; // sample
      }
   groundtruth.close();

   // get thresholds
   ifstream threshold;
   threshold.open("thresholds.txt");
      int th_cloudy, th_xiao_sur_refl_b03;
      double th_msed, th_cos_simi, th_corre_coef;
      threshold >> th_cloudy
         >> th_xiao_sur_refl_b03 
         >> th_msed 
         >> th_cos_simi 
         >> th_corre_coef;
   threshold.close();


   // get range of coordinate to process
   ifstream range2process;
   range2process.open("process_range.txt");
      int top, bottom, left, right;
      range2process >> top >> left >> bottom >> right;
   range2process.close();

   double ndvi_serie_avg[SIZE];
   double ndvi_serie[SIZE];
   int total_cloud;

   // get ndvi serie of pixel is ground truth, calc the avg ndvi serie
   for (int i = 0; i < number_of_gt_pixel; ++i)
   {
      get_px_NDVI_serie(YEAR, gt_coordinate[0][i], gt_coordinate[1][i], th_xiao_sur_refl_b03, ndvi_serie, &total_cloud);
      for (int d = 0; d < SIZE; ++d) ndvi_serie_avg[d] += ndvi_serie[d];
   }
   for (int d = 0; d < SIZE; ++d) ndvi_serie_avg[d] = ndvi_serie_avg[d]*1.0/number_of_gt_pixel;

   cout << "Pre-processing : completed;" << endl;

   int total_cloud_pixel = 0;
   double msed, cs, cc;
   int simi_msed, simi_cos, simi_corr;
   int number_px_padding_rice_msed = 0; 
   int number_px_padding_rice_cs = 0; 
   int number_px_padding_rice_cc = 0;
   ofstream msed_file, cs_file, cc_file;
      msed_file.open("msed.csv");
      cs_file.open("cs.csv");
      cc_file.open("cc.csv");
   // for each pixel in process range, get ndvi, compare to ground truth ndvi
   // also, write re result to file
      // there a log.txt and 3 map in bin text
   for ( int sample = top; sample <= bottom; sample++ )
   {
      for ( int line = left; line <= right; line++ )
      {
         get_px_NDVI_serie(YEAR, line, sample, th_xiao_sur_refl_b03, ndvi_serie, &total_cloud);
         if ( total_cloud > th_cloudy )
         {
            total_cloud_pixel++;
            msed = 100; simi_msed = 0;
            cs = 100; simi_cos = 0;
            cc = 0; simi_corr = 0;
         } else {
            msed = get_Minimum_squared_Euclidean_distance(ndvi_serie, ndvi_serie_avg);
            cs = get_Cosine_similarity(ndvi_serie, ndvi_serie_avg);
            cc = get_Correlation_coefficient(ndvi_serie, ndvi_serie_avg);

            if ( msed <  th_msed ) { simi_msed = 1; number_px_padding_rice_msed++; } else simi_msed = 0;
            if ( cs < th_cos_simi ) { simi_cos = 1; number_px_padding_rice_cs++; } else simi_cos = 0;
            if ( cc > th_corre_coef) { simi_corr = 1 ; number_px_padding_rice_cc++; } else simi_corr = 0;
         }
         cout << line << " " << sample << "\t" << msed << " " << simi_msed << "\t" << cs << " " << simi_cos << "\t" << cc << " " << simi_corr << endl;
         msed_file << simi_msed << ",";
         cs_file << simi_cos << ",";
         cc_file << simi_corr << ",";
      }
      msed_file << endl;
      cs_file << endl;
      cc_file << endl;
   }
   cout << (bottom-top)*(right-left) << endl;
   cout << total_cloud;
   //msed_file << number_px_padding_rice_msed; 
   msed_file.close();
   //cs_file << number_px_padding_rice_cs; 
   cs_file.close();
   //cc_file << number_px_padding_rice_cc; 
   cc_file.close();
   time_t show_end = time(0);
   cout << "All tasks: completed in " << show_end - show_start << " second;\n";
   return 0;
}


void get_px_NDVI_serie(int year, int line, int sample, double th_xiao_sur_refl_b03, double ndvi_serie[], int *total_cloud)
{
   *total_cloud = 0;
   string sql;
   stringstream sqlstream;
      
   try{
      connection C("dbname=mymodisdb user=demo1usr password=demo1pwd hostaddr=127.0.0.1 port=5432");
      if ( !C.is_open() ) { return; }

      sqlstream.str("");
      sqlstream << "\
         SELECT pixels.px_x, pixels.px_y, pixels.date_tkn, pixels.ndvi, pixels.lswi, pixels.evi, pixels.sur_refl_b03 \
         FROM public.pixels \
         WHERE px_x = "<< line <<" and px_y = "<< sample <<" \
         ORDER BY pixels.date_tkn \
         LIMIT 46 OFFSET "<< 46*(year-2009);
      sql = sqlstream.str();
      nontransaction N(C);
      result R( N.exec( sql ));

      int date = 0;
      double ndvi;
      for (result::const_iterator c = R.begin(); c != R.end(); ++c)
      {
         ndvi = 0.0;
         if ( c[6].as<double>() < th_xiao_sur_refl_b03 ){
            ndvi = c[3].as<double>();
         } else {
            (*total_cloud)++;
         }
         ndvi_serie[date] = ndvi;
         date++;
      }
      C.disconnect ();
   }catch (const std::exception &e){
      cerr << e.what() << std::endl;
      return;
   }
}

void clear_array(double a[], int size)
{
   for (int i = 0; i < size ; i++ ){
      a[i] = 0;
   }
}

void print_array(double a[], int size){
   for (int i = 0; i < size; ++i)
   {
      cout << a[i] << " ";
   }
   cout << endl;
}

double get_Minimum_squared_Euclidean_distance( double x[], double y[] )
{
   return get_sum_diff2( x, y);
}

double get_sum_diff2( double x[], double y[])
{
   double total = 0.0;
   for (int i = 0; i <  46; ++i)
   {
      total += ( x[i] - y[i] )*( x[i] - y[i] );
   }
   return total;
}

double get_Cosine_similarity( double x[], double y[] )
{
   return get_sum_product(x, y) / sqrt( get_sum_product(x) * get_sum_product(y) );
}

double get_sum_product(double x[])
{
   get_sum_product(x, x);
}

double get_sum_product(double x[], double y[])
{
   double total = 0.0;
   for (int i = 0; i < 46; ++i)
   {
      total += x[i] * y[i];
   }
   return total;
}

double get_Correlation_coefficient( double x[], double y[] ){
   return get_sum_product_diff_avg(x , y) / sqrt (get_sum_product_diff_avg(x) * get_sum_product_diff_avg(y) );
}

double get_sum_product_diff_avg(double x[])
{
   return get_sum_product_diff_avg(x, x);
}

double get_sum_product_diff_avg(double x[], double y[])
{
   double total = 0.0;
   
   double avg_x = 0.0;
   double avg_y = 0.0;
   for (int i = 0; i < 46; ++i)
   {
      avg_x += x[i];
      avg_y += y[i];
   }
   avg_x = avg_x/46; avg_y = avg_y/46;

   for (int i = 0; i < 46; ++i)
   {
      total += ( x[i] - avg_x ) * ( y[i] - avg_y );
   }
   return total;
}